load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"


fc    = addfile("/mypath/ecmwf/sic_19792004_nov.nc","r")
fw    = addfile("/mypath/ecmwf/sic_20092018_nov.nc", "r")

;*******************************************;
;Read in variables
;*******************************************
lat=fc->latitude
lon=fc->longitude
nlat  = dimsizes(lat)
nlon  = dimsizes(lon)

zc=fc->sic
zw=fw->sic

;***********************************************
; average over the simulations/years
;************************************************

tcs=zc(0,:,:)
tcs=dim_avg(zc(latitude|:,longitude|:,time|:))
tws=zw(0,:,:)
tws=dim_avg(zw(latitude|:,longitude|:,time|:))


;************************
; calculate the anomaly:
;************************

  finfo_ano=tws
  finfo_ano=(tws - tcs)*100.

  finfo_ano!0 = "lat"

  finfo_ano!1 = "lon"
;***********************************
; perform t-test
;***********************************

  xtmp = zc(latitude|:,longitude|:,time|:)     ; reorder but do it only once [temporary]
  ytmp = zw(latitude|:,longitude|:,time|:)

  xAve = dim_avg (xtmp)              ; calculate means at each grid point
  yAve = dim_avg (ytmp)
  xVar = dim_variance (xtmp)         ; calculate variances
  yVar = dim_variance (ytmp)
  sigr = 0.1                        ; critical sig lvl for r
  xEqv = equiv_sample_size (xtmp, sigr,0)
  yEqv = equiv_sample_size (ytmp, sigr,0)
;
   xN=xEqv
   yN=yEqv
   iflag= True                         ; population variance not similar
 ;iflag= False
  prob = ttest(xAve,xVar,xN, yAve,yVar,yN, iflag, False)
  prob!0 = "lat"
  prob!1 = "lon"
  prob&lat =lat
  prob&lon =lon

finfo_ano@_FillValue=-2147483647

do i=0,nlat-1
do j=0,nlon-1
if (.not.ismissing(finfo_ano(i,j)) ) then
if (abs(finfo_ano(i,j)).lt.0.00001)
finfo_ano(i,j)=finfo_ano@_FillValue
end if
end if

end do
end do

;-------------------------
; making the plots
;-------------------------

;-------------------------
wks = gsn_open_wks("ps","NHallanoice_95per_nov_medpolar")

gsn_define_colormap(wks,"BlueDarkOrange18")

plot = new(1,graphic)                          ; create a plot array

  res          = True
  res@gsnDraw  = False                          ; don't draw
  res@gsnFrame = False                          ; don't advance frame, these 2 are going to be set with panel plot - gsn_panel!
  res@cnInfoLabelOn = False                     ; turn off cn info label
  res@cnLinesOn           = False
  res@cnFillOn            = True          ; turn on color
  res@lbLabelBarOn        = True           ; turn off individual cb's

  res@mpPerimOn = True
  res@cnLevelSelectionMode =  "ExplicitLevels"

 res@cnLineLabelsOn=False
 res@mpCenterLonF     = 300.

 res@mpFillOn = True
 res@mpLandFillColor = "grey"
 res@mpFillDrawOrder = "PostDraw"
 res@cnMissingValFillColor   = "white 
 res@gsnPolar   = "NH"
 res@mpMinLatF            = 50.

 res1          = True
 res1@gsnDraw  = False                          ; don't draw
 res1@gsnFrame = False                          ; don't advance frame, these 2 are going to be set with panel plot - gsn_panel!
 res1@cnInfoLabelOn = False                     ; turn off cn info label
 res1@cnLinesOn           = False; True
 res1@cnFillOn            = True          ; turn on color
 res1@lbLabelBarOn        = False           ; turn off individual cb's

res1@cnMonoFillPattern = False ;True
res1@cnFillPatterns      = (/-1,8,-1/) ;17, 8, 13
res1@cnLevelSelectionMode =  "ExplicitLevels"
res1@cnLevels =(/0.0, 0.05/)  ;  90% level

res1@cnMonoFillColor     = True         ; Use same fill color

res1@cnFillColor         = (/"black"/)

res1@cnLineLabelsOn=False
res1@cnMonoLineColor= True
res1@cnLineColor      = (/"grey36"/) ; (/"white"/)

res1@gsnCenterString = ""
res1@gsnLeftString = ""
res1@gsnRightString = ""
res@gsnLeftString = ""
res@gsnCenterString = ""

res@cnFillColors = (/2,3,4,5,7,9,11,12,13,14,15,16/)
res@cnLevels =(/-60.,-40.,-20.,-10.,-5.,-2.5,0,2.5, 5./)

res@gsnRightString = ""
plot(0) = gsn_csm_contour_map(wks,finfo_ano(:,:),res)
plot1 = gsn_csm_contour(wks,gsn_add_cyclic_point(prob),res1)
overlay(plot(0),plot1)

;________________________
; creating the panel:
;________________________
 resP            = True

resP@gsnPanelLabelBar    = False;True                ; add common colorbar
resP@lbLabelFontHeightF  = 0.007               ; make labels smaller

resP@gsnPanelBottom   = 0.05                   ; add space at bottom
;resP@gsnPanelFigureStrings= (/"a)","b)","c)"/) ; add strings to panel
res@txFontHeightF     = .24


resP@gsnMaximize    = True                ; maximize plot
gsn_panel(wks,plot,(/1,1/),resP)
